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ABSTRACT 

Phase-diversity techniques provide a novel observational method for overcomming the 
effects of turbulence and instrument-induced aberrations in ground-based astronomy. Two 
implementations of phase-diversity techniques that differ with regard to noise model, esti- 
mator, optimization algorithm, method of regularization, and treatment of edge effects are 
described. Reconstructions of solar granulation derived by applying these two implementa- 
tions to common data sets are shown to yield nearly identical images. For both implemen- 
tations, reconstructions from phase-diverse speckle data (involving multiple realizations of 
turbulence) are shown to be superior to those derived from conventional phase-diversity data 
(involving a single realization). Phase-diverse speckle reconstructions are shown to achieve 
near diffraction-limited resolution and are validated by internal and external consistency tests, 
including a comparison with a reconstruction using a well-accepted speckle-imaging method. 

Subject headings: methods: observational, techniques: image-processing, Sun: granula- 
tion 
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1. Introduction 

An important goal in ground-based astronomy is to 
improve the angular resolution that can be achieved. 
The angular resolution is nearly always limited by 
phase aberrations introduced by atmospheric turbu- 
lence. For the special case of ground-based solar as- 
tronomy, the spatial resolution is typically limited to 
about O'/S for short-exposure images (< 20 ms) and to 
about l'/0 for long-exposure images (« Is). Many ba- 
sic processes on the sun, however, take place at scales 
below l'.'O. The photon mean free path in the lower 
photosphere corresponds to about O'/I at disk center. 
Magnetic structures may occur on even smaller scales. 
For example, magnetic elements outside of sunspots 
have typical diameters smaller than 0 / /2 (Keller 1995). 
Despite their small size, these small structures are 
believed to play an important role in large-scale phe- 
nomena such as the solar magnetic dynamo or irra- 
diance variations. To understand a variety of solar 
phenomena, it is, therefore, indispensable to reach a 
spatial resolution well below the seeing limit and pos- 
sibly approaching the diffraction limit of existing and 
future, large solar telescopes. 

A number of sophisticated techniques have been 
conceived to combat the deleterious effects of atmo- 
spheric turbulence in astronomical imaging in general. 
Among these are speckle imaging, phase diversity, and 
phase-diverse speckle imaging. In this paper we eval- 
uate the use in solar astronomy of phase-diversity and 
phase-diverse speckle, referred to jointly as phase- 
diversity techniques. To undertake this evaluation, 
we have applied phase-diversity techniques to solar- 
granulation data collected with the Swedish Vacuum 
Solar Telescope (SVST) on La Palma. A subset of 
these data was also processed with a conventional 
speckle-imaging method to demonstrate consistency 
between accepted and novel restoration techniques. 
Phase-diversity techniques are particularly attractive 
for solar astronomy because (I) they require relatively 
simple and inexpensive instrumentation, (2) they per- 
form well with relatively few images in high-signal 
regimes, (3) they lead to a joint estimation of the ob- 
ject and the wavefront, and (4) they obviate the need 
for complicated calibration. 

In the following section we summarize three rele- 
vant fine-resolution imaging techniques: conventional 
speckle imaging, phase diversity, and phase-diverse 
speckle. Phase-diversity techniques, including phase 
diversity and phase-diverse speckle, have been imple- 


mented differently by the Environmental Research In- 
stitute of Michigan (ERIM) group and researchers at 
the Stockholm Observatory that operate the SVST 
(referred to herein as the SVST group). These im- 
plementations are described in Section 3. Section 4 
provides details of the data collection and processing. 
Results derived from applying phase-diversity tech- 
niques to various data subsets are presented in Sec- 
tion 5. A speckle restoration is also included in this 
section to provide an external reference. Conclusions 
regarding the suitability of phase-diversity techniques 
for solar imaging are drawn in Section 6. 

2. Fine-Resolution Imaging Techniques 

Speckle imaging is a relatively mature technique 
for obtaining fine- resolution images in the presence of 
atmospheric turbulence. This technique requires the 
collection of many short-exposure images of a static 
object. The exposure time for each frame must be 
short enough that the evolving atmosphere can be 
regarded as frozen during the exposure. Clever pro- 
cessing of this short-exposure time series affords the 
restoration of fine-resolution information that would 
be irretrievable if a single, long-exposure image were 
collected (Dainty 1984). Speckle imaging requires the 
collection of many images (typically 100 or more) so 
that the ensemble average over the class of all possible 
realizations may be approximated by an arithmetic 
average over a finite number of realizations. Speckle- 
imaging methods have been successfully adapted to 
the solar-imaging problem (Keller & von der Liihe 
1992, de Boer & Kneer 1992, von der Liihe 1993, 
1994). 

Another technique whose aim is the restoration of 
fine-resolution images in the presence of phase aberra- 
tions is the method of phase diversity, first proposed 
by Gonsalves (1979, 1982). Hogbom (1988) indepen- 
dently proposed a special case of this same technique, 
calling it the focal-volume method. Phase diversity 
requires the simultaneous collection of two (or more) 
short-exposure images. Typically, one of these im- 
ages is the conventional focal-plane image that has 
been degraded by the unknown aberrations. A sec- 
ond image of the same object is formed by perturb- 
ing the unknown aberrations in some known fashion 
This can be conveniently accomplished with a sim- 
ple beam splitter and a second detector array that is 
translated along the optical axis. An image collected 
in this second optical channel will contain the effects 
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of the unknown phase aberrations but will also be 
influenced by the intentional defocus, which adds a 
known quadratic phase. It is somewhat remarkable 
that estimates for the object and the unknown aber- 
rations can be made from these two images, given 
the known quadratic phase diversity. The first use 
of the method of phase diversity to retrieve fine- 
resolution solar images was recently reported (Lofdahl 
& Scharmer 1994a, b). 

Intuition suggests that in stressing regimes (eg. 
poor seeing or weak signal levels) a single pair of 
phase-diversity images may not contain enough in- 
formation to estimate jointly and with high fidelity 
the object and wavefront. Even under favorable con- 
ditions for which phase-diversity is able to produce 
good wavefront estimates, we have observed that ob- 
ject information at isolated spatial frequencies may 
be irretrievably lost, resulting in object estimates 
with significant artifacts. These cases motivate a 
third fine-resolution imaging technique, referred to 
as phase-diverse speckle (Paxman, Schulz, & Fienup 
1992a, Paxman & Seldin 1994). As its name sug- 
gests, phase-diverse speckle blends the fundamen- 
tal concepts of phase diversity and speckle imaging. 
Phase-diverse speckle requires the simultaneous col- 
lection of one conventional short-exposure image and 
at least one short-exposure image with phase diver- 
sity, for each of multiple atmospheric realizations, as 
depicted in Figure 1. This makes for a relatively sim- 
ple data-collection scheme. Fortunately, the primary 
strengths of the two constituent methods, namely the 
added information content in a sequence of short- 
exposure images and the wavefront identification pro- 
vided by phase diversity, persist. Two different pro- 
cessing approaches have been demonstrated with real 
phase-diverse speckle data by Lofdahl and Scharmer 
(1994b), referred to here as partitioned phase-diverse 
speckle (PPDS, see section 3.2.) and by Seldin and 
Paxman (1994), called joint phase-diverse speckle 
(JPDS, see section 3.1.). 

In order to draw precise distinctions between these 
fine-resolution imaging methods and to establish no- 
tation, we present our working data-collection model. 
We concentrate on the data-collection model for phase- 
diverse speckle, from which it can be seen that phase- 
diversity and speckle-imaging data sets are special 
cases. The incoherent isoplanatic image-formation 
process is well modeled by the following discrete con- 
volution: 


9jk{x) = ^Tf(x')s jk {x-x'), (1) 

X> 

where f(x) is the object array, sjk(x) is an incoherent 
point spread function (PSF) corresponding to the j th 
atmospheric realization and the Arth diversity channel, 
gjk(x) is the corresponding noiseless image, and x is a 
two-dimensional coordinate. The size of the noiseless 
image array is determined by the field of view (FOV) 
of the detector, whereas the size of the object array 
should extend well beyond the FOV. Of course any 
detected images will contain noise. The detected data 
set is represented by {djk}, where 

i : I’ 2’ ■ ( ' 2) 

and where the general noise operator Af[ *] introduces 
photon noise, additive Gaussian noise, and/or any 
other noise sources that are appropriate. The data 
set contains image frames from a total of J atmo- 
spheric realizations and K diversity channels, where 
typically K = 2. Phase diversity is introduced by 
including a known phase function in the system's co- 
herent transfer function (Goodman 1968). 

Pjk(u) = P{u)exp{i[<f)j{u) + 6k(u)}} , (3) 

where P(u) is a binary function that serves as an ap- 
propriately scaled model of the telescope pupil func- 
tion, ) is an unknown phase-aberration function 
with contributions from the jth atmospheric realiza- 
tion and the fixed telescope aberrations, 0*(u) is a 
known phase-diversity function associated with the 
kth diversity channel, and u is the discrete spatial- 
frequency variable. The phase-diversity function. 
0*(u), will be zero in the conventional channel and 
quadratic in the channel with defocus. It is conve- 
nient to parameterize the phase-aberration function 
using coefficients for an appropriate set of basis func- 
tions (such as Zernike polynomials): 

Af 

= ^ ^ * ( 1 ) 

m=l 

The incoherent PSF, Sj*(.r), is just the squared mod- 
ulus of the discrete Fourier transform (DFT) of the 
coherent transfer function in equation (3). Thus the 
noiseless image, gjk{x), implicitly depends upon both 
the object and aberration parameters. 

The goal of phase-diverse speckle is to estimate t he 
common object and each of the J phase-aberration 
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functions (or equivalently the aberration parameters) 
from the J K detected images, given the known phase- 
diversity functions and the binary pupil function. 
Notice that this model accommodates conventional 
phase diversity when J — 1. When K = 1, the 
dataset corresponds to conventional speckle data. Al- 
though conventional speckle processing seeks an ob- 
ject estimate, no attempt is made to estimate the in- 
dividual phase-aberration functions associated with 
the atmospheric realizations. 

3. Implementation of Phase-Diversity Tech- 
niques 


eters for each atmospheric realization, is given by 


K 


■mw*}] = n n n 

j=l k=l x r 


( 5 ) 

We use the principle of maximum likelihood to solve 
the inverse problem. We jointly estimate the object 
and aberration parameters by maximizing the log of 
the likelihood function, 


j K 

£(/.«) = • * ,5) 

j = 1 k = 1 x 


Both the ERIM and SVST research groups have 
been working on phase-diversity techniques for sev- 
eral years. Although the data-collection paradigm 
and basic goals are common to both groups, the pro- 
cessing implementations differ, reflecting differing re- 
search paths, emphases, and insights. In this section 
we summarize the salient features of these differing 
implementations and quote references that provide 
details of the processing. 


where a constant term that has no bearing on the 
maximization procedure has been dropped and the 
phase-aberration parameter estimates, a jm , have been 
lexicographically arranged into a single vector, a. We 
use the caret symbol,^, throughout to indicate an es- 
timate. Because aberration parameters for all J real- 
izations are estimated simultaneously along with the 
object parameters, we refer to reconstructions as joint, 
phase-diverse speckle (JPDS) estimates. 


3.1. ERIM Implementation 

A guiding philosophy of the ERIM group has been 
to model the forward problem (data collection) as ac- 
curately as possible and to use this model as the basis 
for solving the inverse problem (object and aberration 
estimation) using estimation-theoretic tools. 

3.1.1. ERIM noise model and likelihood function 

A Poisson noise model was selected because such 
a model accurately accommodates the combined ef- 
fects of signal-dependent photon noise and additive 
Gaussian CCD readout noise. The number of photo- 
conversions that occur at each detector element will 
be a Poisson-distributed random variable with a mean 
value prescribed by the noiseless image, gjk(x) y given 
in units of mean detected photons per pixel. Although 
not explicitly shown here, an artificial bias is added to 
the noiseless image to model the readout noise (Sny- 
der, Hammoud, & White 1993). Assuming that the 
detected signal at each detector element is statisti- 
cally independent, the probability of acquiring a data 
set {djk}, given the object and the aberration param- 


3.1.2. ERIM optimization algorithm 

The method of preconditioned conjugate gradi- 
ents (Luenberger 1984), a conventional nonlinear- 
optimization technique, is employed to maximize equa- 
tion (6) over the set of object pixel values and phase- 
parameters. Conjugate-gradients optimization is an 
iterative procedure that, at each iteration, requires a 
single gradient computation and a line search involv- 
ing repeated likelihood evaluations. A closed- form 
expression for the gradient of the log-likelihood func- 
tion has been derived (Paxman, Schulz, & Fienup 
1992b, Paxman & Seldin 1994) and is used extensively 
in the iterative search. 

3.1.3. ERIM regularization 

The maximum-likelihood estimates of the object 
pixels can be somewhat sensitive to noise and may 
require a regularization strategy in which resolution 
in the estimate is traded for noise suppression. There 
are many candidate regularization strategies, however 
in this implementation the method of sieves (Snyder 
& Miller 1985) is employed. This is accomplished by 
constraining the object estimate to be of the form 

/(*) = /(*')«(* - x ') ’ < 7) 

X 1 
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which is the convolution of an artificial Poisson point 
process, f(x ), with a smoothing (or sieve) kernel, 
t?(or). The maximum-likelihood formulation remains 
unchanged. However, instead of estimating the ob- 
ject, we now estimate the new underlying process, 
f(x). The final object estimate is formed from f(x) 
and t»(x) using equation (7). The choice of an ap- 
propriate sieve is the subject of ongoing research, but 
the Gaussian kernel (Snyder fe Miller 1985) has been 
quite effective. 

3.1.4 . ERIM treatment of edge effects 

In principle, the size of the FOV is limited by the 
extent of the detector array or a field stop. How- 
ever, the N x N FOV over which the convolutional 
imaging model in equation (1) is valid will depend on 
anisoplanatic effects. We therefore treat an N x N 
subframe of data as an effective FOV (as if it were 
collected by a detector array of that extent), and we 
estimate the object-pixel values associated with these 
elements. In addition, we estimate object-pixel values 
within a guard band of width B pixels surrounding 
the effective FOV. Although these guard-band pixels 
do not have a corresponding detector element within 
the effective FOV, they influence the data in two dis- 
tinct ways. The first is through PSF sidelobes. For 
example, a bright object point in the guard band will 
create PSF sidelobes that spill into the effective field 
of view. The second mechanism derives from random 
image translations that occur as a result of differing 
tilt components among the atmospheric realizations 
in a phase-diverse speckle data set. Thus the main 
lobe of a PSF associated with an object point in the 
guard band may be directly sensed by detector ele- 
ments within the effective FOV when the tilt com- 
ponent for a particular realization provides the right 
translation. The size of the guard band is selected so 
that the influence of pixels far from the effective FOV 
is negligible. This is related to the severity of the 
aberrations and the resulting PSF side-lobe structure 
and translations. The total number of object param- 
eters is given by (N 4* 2 B) 2 . 

Several aspects of the guard-band technique are 
appealing. The guard-band technique affords the ac- 
curate and efficient computation of the convolution 
in equation (1) using a DFT, or a fast Fourier trans- 
form (FFT) in practice. Although a DFT assumes 
a periodic object, estimated values for pixels at the 
outer rim of the guard band are allowed to “wrap 
around” since their influence on the estimated data 


is negligible. Another appealing aspect of the guard- 
band method is that, unlike apodization techniques 
(Paxman h Crippen 1990), the measured data are 
unperturbed. Finally, the guard-band method allows 
for the reliable retrieval of object pixels up to the edge 
of the detector-limited FOV, so long as the effective 
FOV is defined to abut the detector-limited FOV. 

3.2. SVST Implementation 

The purpose of the SVST-group implementation 
is to develop a fast and reliable method for obtain- 
ing nearly diffraction-limited images with the SVST 
in La Palma. The code has been operational since 
the spring of 1993 and is the first phase-diversity 
code used to demonstrate, through several consis- 
tency tests (Lofdahl & Scharmer 1994a, b), that the 
technique works on real data. This is mainly due to 
the successful implementation of methods to deal with 
image boundary effects and for registration of focused 
and defocused images pairs. 

3.2.1. SVST noise models and metric 

The SVST code is based on two modifications of 
the error metric of Gonsalves and Chidlaw (1979). 
which relies on a Gaussian additive noise model. This 
model allows the estimation of the optimum object to 
be performed implicitly while the estimation of the 
optimum wavefront is done explicitly, which leads to 
a straightforward and computationally efficient code. 
However, the expression for the optimum object can- 
not be used directly because it leads to unlimited am- 
plification of noise at spatial frequencies where the 
transfer functions of the focused and defocused im- 
ages simultaneously approach very small values. This 
happens because the expression for the restored ob- 
ject gives a best fit to the data, including its noise. 
Of course what is needed is an expression for the re- 
stored object which is as accurate as possible, which 
means that the observed images must be filtered to re- 
duce noise in the restored object. We have also found 
that high noise levels give slower convergence in the 
iterative determination of the wavefronts. 

Our first modification to the error metric of Gon- 
salves and Chidlaw, therefore, is to introduce a noise 
filter, which is applied both to the observed focused 
and defocused images, in the expression for the error 
metric. 

The second modification consists of accounting for 
the possibility that the signal-to-noise ratio (SNR) 
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may be significantly different for the image channels, 
as in the case of beam-splitters that do not distribute 
light in equal proportions. With these two modifica- 
tions the error metric becomes 


Lj = ^\H j D jl -F j S jl \ 2 +-i\H j D j2 -F j S j2 \ 2 , (8) 

li 

where, in contrast to the ERIM formulation, the sum- 
mation is in the Fourier domain rather than the im- 
age domain. Sjk is an estimate of the optical transfer 
function (OTF), which is the Fourier transform of the 
estimated PSF. Hj is the noise filter, and 7 is given 
by 

7 = v\l ° 2 . (9) 

where 07 and cr? are the RMS values of the noise of 
the two channels. With Hj = 1 and 7=1, the error 
metric of Gonsalves and Chidlaw (1979) is recovered. 
Following their derivation, which means estimating 
F independently for each realization j, leads to an 
expression for the estimated Fourier object, 


f. — ft Ril^h + 2 
J ~ j | 5 jl | 2 + 7 | 5 j2 | 2 


( 10 ) 


Hj , but refers now to several atmospheric realizations. 
These filters are specified in Section 3.2.3. 

Because phase estimates derive from partitioned 
data whereas object estimates derive from these phase 
estimates in conjunction with a full phase-diverse 
speckle data set, we refer to this processing approach 
as partitioned phase-diverse speckle (PPDS). 


3.2.2. SVST optimization algorithm 

The expansion for (f>j , equation (4), allows us to 
write Sjk = Sjk{ocj m ) y and therefore Ej = Ej(a jm ). 
Due to the nonlinear dependence of Ej on a jm . the 
minimum has to be found iteratively from an initial 
estimate, usually <f>j = 0. This is implemented by 
approximating changes in Ej by 




(H) 


and seeking corrections to the coefficients, Saj m , such 
that the minimum of Lj is found in the next iteration. 
This linearization results in a matrix equation of the 
type 

A • + b = 0 , ( 1 0 ) 


where * used as a superscript denotes the complex 
conjugate. The corresponding error metric can be 
written as 

£j=E \Ej\ 2 , (11) 

u 

where the Fourier-domain error function is defined as 

Ej = Hj d ; 2S}1 ~ Mr , ( 12 ) 

Vl%l 2 + 7|5,- 2 | 2 


where the elements of A and b are sums of different 
combinations of Ej and its partial derivatives with 
respect to the aberration parameters (see Lofdahl fe 
Scharmer 1994b). These derivatives involve the trans- 
formed images, Djk , the OTFs, Sjk , and the OTF 
derivatives, which are evaluated analytically. Note 
that A is an M x M matrix, where M is the num- 
ber of aberration parameters. This equation is solved 
with the singular value decomposition (SVD) method, 
as described in Section 3.2.3. 


Earlier analysis (Lofdahl & Scharmer 1994b) has 
shown that this method leads to good estimates of 
the wavefront parameters but that in poor seeing the 
restored objects are contaminated by artifacts from 
poor SNR at isolated spatial frequencies. These ar- 
tifacts are removed by combining the results of two 
or more realizations to calculate the restored Fourier 
object in a fashion that is well known in the literature, 


F = 


H 


DjiSU 

£/=ilS,-i| 2 +7|5 i2 | 2 


(13) 


3.2.3. SVST regularization 

In the SVST implementation, two methods of reg- 
ularization (noise suppression) are employed: (l) the 
observed images are low-pass filtered and (2) the 
wavefront estimate is restricted by zeroing the least 
significant singular values of the system matrix, A. 

In order to provide noise reduction, the filters // } 
and H must be specified. Since the main priority is 
to obtain good, restored images, it seems reasonable 
to choose H such that the combined RMS error from 
noise and the filter is minimized in the restored objec t , 


where the OTFs must include the shifts necessary to 
bring all images into co-alignment. The filter H is of 
similar form to the filters of individual realizations, 
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F. This leads to 


H = l-{\N 1 \ 2 } 


E/ =1 l$il a + 7 & 2| 2 \ 

I Ej=i DjiSji + 7 HjoS* o J 


( 16 ) 

where {■) denotes an expectation value. In prac- 
tice, the second expectation value is estimated with a 
smoothing operation and the removal of noise peaks in 
the high frequency regime, see (Lofdahl & Scharmer 
1994b) for details. 

The filter area expands automatically with the 
number of included realizations, as information is 
added in F at isolated frequencies and as the noise 
influence is reduced at high frequencies by the av- 
eraging of many realizations. The single-realization 
filters, Hj , are defined as the special case where the 
sum is over one realization. 


Solving the matrix equation (15) directly for a 
large number of wavefront parameters gives poor con- 
vergence, or even divergence, in the iterative proce- 
dure and poor wavefront estimates. One conjecture 
is that this happens because the focused and defo- 
cused images do not contain enough information to 
distinguish between too many wavefront parameters. 

The matrix equation is therefore solved by means of 
the SVD method. The SVD method rearranges the 
equations, so that solutions are sought for orthogonal 
linear combinations of the parameters. These com- 
binations are sorted in order of significance, as ex- 
pressed by the singular values. Solving only the sys- 
tem of equations with singular values larger than a 
cutoff level, defined as a fraction of the largest sin- 
gular value, restricts the solution to the subspace 
spanned by the most significant linear combinations, 
corresponding to the retained equations. 

Recent experiments within the SVST group have 
shown that using a cutoff level of 0.02 permits Zernike 
parameters up through the 12th radial degrees and 
azimuthal frequencies and generates wavefronts that 
conform better to Kolmogoroff covariance. This value 
of the cutoff level was found by trial and error to give 
only a 1% increase in the converged value for the error 
metric L . The SVD method eliminates the problem 
of over-parameterization in a simple and automatic 
way, independent of the type of object used to deter- 
mine the wavefront parameters. It also significantly 
improves the convergence of the iterations, thus en- 
hancing the speed of the code. 

With the low-order wavefront expansion used for 
the current work, the inversion problem is w'ell-conditioned. 


We therefore used a cutoff of 0.0001, which effec- 
tively de-activates the regularizing property of the 
SVD method. 

3,2.4- SVST treatment of edge effects 

The error metric of Gonsalves, as expressed in the 
Fourier domain, does not include the effects of bound- 
aries of the images. Using FFTs to perform the con- 
volutions needed for the calculation of the error func- 
tions and its derivatives produces severe wrap-around 
effects which can lead to very large errors in the de- 
rived wavefronts when using small sub-fields. This 
problem is avoided by transforming the error function 
and its derivatives to the image plane. ParsevaFs re- 
lation permits the summation in equation (11) to be 
calculated in the image domain. Observing that the 
wrap-around errors in the image-domain error func- 
tion, ej, are concentrated to the boundaries of the 
array, we use an array size that is large enough that 
the wrap-around effects are accommodated outside 
the effective FOV. The summation is then restricted 
to the FOV. Apodizing with a modified Hanning win- 
dow function with a flat profile over the FOV further 
removes a high frequency pattern from the disconti- 
nuities at the array boundaries (Lofdahl fe Scharmer 
1994b). 

Like the ERIM algorithm, the technique affords 
the accurate and efficient computation of the con- 
volution in equation (1) using FFTs. Furthermore, 
unlike pure apodization techniques (Paxman Crip- 
pen 1990), the measured data are unperturbed in the 
N x N area. 

4. Solar Data 
4.1. Data Collection 

The data were collected with the 47.5 cm SVST in 
La Palma on April 27, 1993 (see Lofdahl & Scharmer 
1994b for details). The multi-image, real-time frame 
selection system (see Scharmer & Lofdahl 1991) mon- 
itored the fine-structure content in the focused chan- 
nel and was set to select and store the best 100 image 
pairs out of 1500 recorded during 30 second intervals. 
The best seeing at the SVST usually occurs intermit- 
tently, so the series used for this analysis was recorded 
during 6 seconds of good seeing at 14:19 UT. This in- 
terval is short enough that no significant evolution of 
the granulation structure can take place. Observa- 
tions were made through a 5.4 nm wide interference 
filter centered at 470 nm. The known quadratic phase 
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difference of the two image channels corresponds to 
a phase shift at the edge of the aperture equal to 
0.985 ± 0.06 waves. 

Images were recorded by two synchronized EEV 
video CCD cameras operating at 50 Hz (20 ms ex- 
posure time) and digitized to 8 bits by two Kontron 
DEC/IPS image processing systems. The image scale 
is 0'/ 0706 and 0 / /0732 per pixel in the x and y direc- 
tions respectively. (The mean value, 0'/0719 was used 
for the inversions.) The image scales of the two image 
channels differ by less than 0.1%, and the images are 
rotated by less than 0.1 degrees relative to each other. 


4.2. Preprocessing 

The frames in the two channels were corrected with 
their corresponding dark current and gain table cali- 
brations. The gain table was determined by randomly 
moving the telescope and averaging a large number of 
frames. Since the bias level of the cameras changed 
with illumination, the bias was determined from cov- 
ered parts of the CCD sensor and taken into account 
in the data reduction. Image-restoration techniques 
can be sensitive to small but consistent errors in the 
gain table. Despite the excellent results reported by 
Lofdahl and Scharmer (1994b), a thorough analysis 
of the sensitivity of phase-diversity restorations to 
camera calibration errors has not been undertaken. 
Therefore, a careful calibration procedure was fol- 
lowed here to mitigate the influence of calibration 
errors on our analysis. Following the procedure in 
(Keller, Stenflo, & von der Liihe 1992), the gain table 
correction was improved by decomposing the Fourier 
transform of the average image into a high-frequency 
domain and a low-frequency domain. It is then as- 
sumed that the signal in the high-frequency domain 
is due to errors in the gain table only. The true av- 
erage image may then be found with an appropriate 
low-pass filter. Each frame in the sequence is then 
corrected by multiplication with the ratio between 
the low-pass filtered and the original average frame 
according to 


i corrected 
a jk 



(17) 


where ~ represents the spatial low-pass filtering. The 
focused and defocused channels are treated sepa- 
rately. 

Each defocused image was further adjusted so that 
its mean intensity was the same as its corresponding, 


focused image. A close inspection of the defocused 
images revealed a very weak, horizontal strip-pattern 
that could not be removed with the gain-table calibra- 
tion. This pattern was removed in the following way: 
for each image the average column was determined 
by averaging along the horizontal direction. A high- 
order Legendre polynomial was subtracted from the 
average column to extract the high-frequency strip 
pattern and to remove any variations due to solar 
structures. This difference was then subtracted from 
every column in an image. 

Shifts between consecutive frame pairs were deter- 
mined via cross-correlation and then removed. The 
same shifts were applied to both channels, and only 
full pixel shifts were performed in order not to af- 
fect any high-spatial-frequency solar signal. Note that 
this prealignment was performed to accommodate the 
speckle reconstruction and is not required for phase- 
diverse speckle. Finally, 128 x 128 pixel subframes 
were extracted in such a way that the average shift 
between focused and defocused subframes amounts to 
less than a pixel. 

The entire sequence of 100 preprocessed image 
pairs is publically available and can be obtained by 
contacting the authors at the Stockholm Observatory 

4.3. Processing 

The 128 x 128 subframes were selected with a cen- 
ter region of 70 x 70 pixels that contains fine image 
details that are useful for assessing restoration fidelity. 
This central region, or effective FOV, corresponds to 
a 5'/0 x 5'.'0 patch, which is small enough to satisfy an 
isoplanatic imaging assumption (Lofdahl & Scharmer 
1994a, b). Both the conventional speckle and SYST 
group restorations use all the data in the larger 128 
x 128 subframes, but in both cases the restorations 
are most reliable over the center region due to edge 
effects. The ERIM restoration procedure does not use 
data outside the 70 x 70 center region, but the larger 
subframes were utilized in some cases to obtain a bet- 
ter initial object estimate within the guard band of 
size B = 21 pixels. All comparisons between restora- 
tions are made over this central area, and henceforth 
we will restrict our attention to this region. 

Five examples out of the sequence of 100 pairs of 
preprocessed subframes are found in Figure 2. ! he 

top row contains conventional, focused images m i he 
5'/0 x 5'.'0 central region, and the defocused counter- 
parts are found below. The RMS contrast (defined 
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here as the ratio of the standard deviation to the mean 
intensity) averaged over the entire sequence is 7.6% 
and 4.8% for the focused and defocused images, re- 
spectively. Examples of very good seeing are shown in 
Figures 2(a) and (b). Figure 2(c) is an average case, 
and (d) and (e) are examples of aberrated images that 
yield poor conventional phase-diversity restorations. 

It is interesting to note the dramatic effect of defocus 
on the image in (d). 

To facilitate comparisons of restorations, both the 
ERIM and SVST groups parameterized the phase 
aberrations with Zernike polynomials 4 through 15. 
From previous experience we knew that this level of 
wavefront parameterization leads to a well-conditioned 
inversion problem. The two tilt components, Zernike 
polynomials 2 and 3, provide sub-pixel positioning of 
the estimated PSFs. For the case of JPDS, these pa- 
rameters are included in the joint-estimation process. 
For PPDS, they are derived in post- processing by 
cross-correlating the phase-diversity object estimates. 

If the data were not prealigned, then the use of these 
polynomials would be even more critical. Another 
important parameter is the registration between the 
conventional and diversity channels. Typically, the 
ERIM approach is to estimate these parameters along 
with the aberrations and the object. However, to keep 
the implementations as similar as possible we decided 
to use predetermined values. The technique used here 
for estimating these parameters from the average of 
the data in the two imaging channels is discussed by 
Lofdahl and Scharmer (1994b). Other system param- 
eters, such as the amount and the sign of diversity 
and the pixel spacing, were the same for both imple- 
mentations as well. 

5. Results 

In this section we compare estimates made with: 
(1) ERIM phase-diversity techniques, (2) SVST-group 
phase-diversity techniques, and (3) conventional speckle 
imaging. We show that the two implementations of 
phase diversity yield virtually identical object and 
aberration estimates, but that these object estimates 
are susceptible to artifacts, indicating a clear need 
for phase-diverse speckle estimation. We show quan- 
titatively that JPDS is slightly better than PPDS. 
The internal consistency of the phase-diverse speckle 
estimates is demonstrated by analyzing estimates de- 
rived from disjoint data sets, and external validation 
is obtained by comparing these restorations with the 


speckle restoration. Finally, the influence of anisopla- 
natism is considered. 

5.1. Implementation Invariance for Phase Di- 
versity 

The ERIM and SVST group’s approaches to phase- 
diversity techniques are quite distinct. The noise 
model, estimator, optimization algorithm, regulariza- 
tion, and edge treatment are some of the ways in 
which the techniques differ. Despite these differences, 
we have discovered strong evidence to indicate that in 
most cases the different implementations yield, for ail 
practical purposes, identical solutions. The first evi- 
dence of this is presented here for the case of conven- 
tional phase-diversity; i.e., restoration from a single 
pair of images. 

Figure 3 contains the conventional phase-diversity 
object estimates derived from each of the 5 pairs of 
images in Figure 2. The data in Figure 2 were selected 
because the restorations derived from them offer a 
nice comparison of the two implementations. The first, 
trend to note in Figure 3 is that the ERIM restora- 
tions in the top row have slightly higher spatial fre- 
quency content than the SVST group s restorations 
in the bottom row, which is a direct consequence of 
the different approaches to object regularization. T lie 
attempt to restore more fine detail can have the un- 
desirable side-effect of boosting high-frequency arti- 
facts, creating a slightly mottled appearance. Typi- 
cally, when the SVST group’s noise filter is applied 
to an ERIM restoration, the resulting estimate is vi- 
sually indistinguishable from the SVST restoration 
Thus, the overall features in the restorations from the 
two implementations are often quite similar, but the 
SVST group’s restorations appear to be low-pass fil- 
tered versions of the ERIM estimates. 

A useful measure of similarity between two esti- 
mates, fi(x) and / 2 (x), is via the normalized RMS 
error, e, defined via 

where x r brings / 2 into registration with f \ , and av- 
eraging is done over N 2 pixels. The misregistration. 
x r) is estimated to sub-pixel accuracy by interpolat- 
ing the peak of the cross-correlation of the two esti- 
mates. Phase-diversity estimates were made for the 
entire sequence of 100 image pairs, and the average 
error between the two implementations is t 7 = 1.5%. 
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with a standard deviation of 1.0%. Figures 3(a) - (c) 
contain estimates for which the match is very good, 
with errors of 0.8%, 1.2%, and 1.2%, respectively. 
Figure 3(a) is an example of a high-quality phase- 
diversity restoration because it compares favorably 
with phase-diverse speckle restorations presented in 
subsequent sections (see Figure 8). Visual cues can 
be taken from the central portion of the restorations, 
where there is a narrow intergranular lane and a small, 
bright feature at the tip of the arrow. This small fea- 
ture is slightly brighter in the ERIM restoration, re- 
flecting a higher concentration of energy due to the 
recovery of higher spatial frequencies. On the other 
hand, mottle can be observed on some of the larger 
granules. When the SVST-group noise filter for each 
of the restorations in Figures 3(a) - (c) is applied to 
the corresponding ERIM restoration, c drops to 0.4%, 
0.9%, and 1.0%, respectively, the mottle disappears, 
and the differences in the restorations become visually 
imperceptible. 

The restorations in Figures 3(b) and (c) degrade 
somewhat with respect to (a), but what is also no- 
table is that the ERIM and SVST group’s restorations 
share the same features, even when they are false. For 
example, the small features near the narrow inter- 
granular lane at the center of the image are smeared 
into nearby granules in both of the restorations in 
(c). Referring back to Figure 2(c), it is clear why 
these features, which are blurred severely in the orig- 
inal data, are not fully recovered. Careful inspection 
of Figures 3(a) - (c) reveals that, despite differences 
in spatial frequency content, features such as granule 
shape and intensity variations along the granule edges 
are virtually identical. 

Greater differences between implementations are 
observed in Figures 3(d) and (e), with c = 4.0% and 
2.1%, respectively. The SVST group restoration is 
superior to the ERIM restoration in (d), particularly 
with respect to the large granule in the lower left. 
Conversely, the ERIM restoration in (e) does not dis- 
play the stripe artifact that runs at an angle through 
the SVST restoration. Neither of these restorations 
is particularly good, however. The trend observed 
across the 100 restorations is that large discrepan- 
cies between implementations are observed rarely and 
only for cases when the blurring is too severe for con- 
ventional phase diversity to be effective. 

The other major component of phase-diversity restora- 
tion is the estimation of the phase aberration. In 
Figure 4 we show 12 scatter plots of Zernike coef- 


ficient estimates from the two implementations for 
polynomials 4 through 15. In each sub-graph the co- 
efficients estimated by the two implementations for 
each of the 100 image-pairs are scattered about the 
line y = x. We note immediately the high correla- 
tion between the aberrations derived from the two 
implementations. Aside from one or two outliers, the 
aberration estimates are very consistent. One mea- 
sure of agreement is formed by taking the square root 
of the mean over the pupil of the average, squared 
wavefront difference (averaged over 100 realizations). 
This measure of RMS difference between implemen- 
tations is a negligible 0.043 wave. This is well below 
the well-known Marechal aberration-tolerance condi- 
tion of 1 /14th wave RMS phase error (Born & Wolf 
1980). Systems that meet the Marechal condition are 
considered to be well-corrected, producing imagery 
for which the degradation would be difficult to per- 
ceive, Thus, the aberration estimates from the two 
implementations yield point-spread functions that are 
visually indistinguishable. 

5.2. Partitioned vs. Joint Estimation in Phase- 
Diverse Speckle 

Aside from implementation invariance, important 
conclusions to draw from the results presented in Fig- 
ure 3 are that even an above-average conventional 
phase-diversity restoration like the one in (c) still suf- 
fers from residual blur, and that in some cases the 
restorations are dominated by artifacts. This obser- 
vation was made by Lofdahl and Scharmer (1994b) 
and led to a partitioned phase-diverse speckle (PPDS) 
estimation strategy (Section 3.2.1.) in which pairs of 
restorations were combined to produce much higher 
quality object estimates. Similarly, Seldin and Pax- 
man (1994) demonstrated that restorations from the 
same data using a joint phase-diverse speckle (JPDS) 
algorithm (Section 3.1.1.) became better and more 
consistent as image pairs were added. From these 
results it is clear that there is a distinct advantage 
to using phase-diverse speckle instead of conventional 
phase diversity. Given the need for phase-diverse 
speckle, a natural question to address is whether there 
is an advantage to using the joint-estimation approach 
as opposed to a partitioned technique. 

As with conventional phase diversity, we find that 
object estimates derived from the JPDS and PPDS 
algorithms are virtually identical, with a negligible 
difference of e = 0.3% when the entire sequence of 100 
image pairs is used. Not surprisingly, the JPDS and 
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PPDS restorations are also visually identical, as ob- 
served in Figures 8(b) and (c), respectively. Intuition 
suggests that JPDS would provide better estimates 
of both the object and the aberrations because the 
ratio of the number of measurements to parameters 
increases significantly with the addition of each new 
pair of aberrated images of the same object. How- 
ever, the JPDS and PPDS restorations are so similar 
for these data that we are left to wonder if there is a 
distinct advantage to a joint-estimation approach in 
this regime. 

The object estimates alone do not tell the whole 
story, and so the accuracy of the phase aberration 
estimates must also be considered. Given two aber- 
ration estimates, one from conventional phase diver- 
sity and the other from JPDS, we seek to quantify 
the accuracy of each. To do this without the aid of 
a simulation, we generate an estimated point-spread 
function, sjk, from each wavefront, convolve it with 
the corresponding object estimate, and compare the 
resulting image estimate to the measured image. We 
define a fidelity metric for the j th wavefront estimate, 
<f)j , in a normalized mean-squared error sense as 

R 7 ($j) = ZM i fUl (19) 

[i^Ex/Wl 2 

where averaging is done over N 2 pixels, and gjk is 
formed from / and the estimated point-spread func- 
tion, Sjk , which is a function of <$>j, via the convolution 

&*(*) = • ( 20 ) 

X 1 

Note that the numerator of equation (19) is iden- 
tical (aside from a sign change) to the log-likelihood 
function for phase diversity under an additive Gaus- 
sian noise model (Paxman, Schulz, & Fienup 1992b). 
Smaller values of R imply a better match of the 
object and aberration estimates to the data. We 
evaluate R{<t>j) for j = 1,2, — -,100 for both the 
ERIM JPDS wavefronts estimates and the SVST- 
group phase-diversity wavefront estimates. When 
evaluating R for the JPDS aberration estimates, the 
100-realization JPDS restoration in Figure 8(b) is 
used for /. Similarly, the SVST-group conventional 
phase-diversity aberrations estimates are evaluated 
using the 100-realization PPDS object estimate in 
Figure 8(c). 

Figure 5 is a scatter plot of R(<f>j) for PPDS ver- 
sus R(<t>j) for JPDS. The mean fidelity metric for 


the JPDS wavefronts is 1.3%, compared with 1.5% 
for phase diversity. Since the object estimates used 
to generate the fidelity metric differ by only 0.3%', 
we can conclude that most of the error is due to 
sources other than the object. In 97 of 100 cases, 
the fidelity metric for the JPDS aberration estimate 
is less (better) than the wavefront derived with con- 
ventional phase diversity, and the fidelity metrics are 
very close for the other 3 cases. There is a fairly 
strong linear correlation that indicates that the rise 
and fall in fidelity metric tracks fairly well for the 
two cases. Nonetheless, it is clear that there are 
15 to 20 cases for which the phase-diversity fidelity 
metric is significantly worse in a relative sense. In 
an absolute sense we are observing only very small 
improvements in the aberration estimates, which im- 
proves the fidelity metric by only a fraction of a per- 
cent. Apparently, phase-diversity wavefront estimates 
are so good in this regime that there is little room 
for improvement when using JPDS. Even so, it re- 
mains necessary to use multiple realizations to ob- 
tain a good object estimate via PPDS. The reason 
for this is simply that the OTF for any one realiza- 
tion tends to reach very low values at isolated spatial 
frequencies which vary from realization to realization. 
This explanation is consistent with the ringing arti- 
facts found in Figures 3(d) and (e). The fact that, the 
JPDS wavefront estimates are almost always better 
than the PPDS counterparts shows that the joint- 
estimation approach is a successful concept which is 
likely to be superior to PPDS in different imaging 
regimes; e.g. reduced signal levels, smaller isopla- 
natic patches, stronger aberrations, more aberration 
parameters, etc. 

5.3. Internal Consistency of Phase-Diverse 
Speckle Estimates 

Both conventional phase diversity and phase-diverse 
speckle imaging have been investigated extensively 
with past simulation studies. These simulations are 
an important component in establishing the credi- 
bility of phase-diversity image restoration. An even 
more important step was taken by using real solar 
granulation data to demonstrate the consistency of 
aberration estimates from neighboring image patches 
(Lofdahl & Scharmer 1994b). A measure of inter- 
nal consistency with respect to object estimates was 
demonstrated with the same data set for the case of 
JPDS (Seldin k, Paxman 1994) applied to a small 
number of realizations. In this case several JPDS 


11 



restorations were made using different sets of im- 
ages collected within seconds of each other, and these 
restorations were shown to be highly consistent with 
each other. We explore this in more detail here using 
many more realizations, to further demonstrate the 
consistency of object estimates. 

Object estimates are a random process with a mean 
and a variance, and any one restoration is a sample 
drawn from such a process. Like other random pro- 
cesses, the mean of an object estimate can be esti- 
mated using a sample mean derived from independent 
trials. The same is, of course, true for estimating the 
variance. We define our trials as restorations of the 
same object with different input data. For example, 
we formed a conventional phase-diversity restoration 
for each of the 100 image pairs, and from these 100 
restorations we computed a sample mean and vari- 
ance. We also formed samples of 2-realization JPDS 
restorations by partitioning the sequence of 100 pairs 
into 50 disjoint sets, containing J — 2 image pairs 
each. The same partition of the sequence into dis- 
joint sets was done for the cases of J = 5, 10, 25, 
and 50 realizations, and restorations were made for 
every set. Examples of one restoration for each of 
the 6 cases investigated are shown in Figure 6. The 
1- and 2-realization restorations in Figures 6(a) and 
(b), respectively, have visible artifacts, but restora- 
tion quality is quite good in the 5-realization case in 
Figure 6(c). It is difficult to see much variation across 
the 10-, 25-, and 50-realization restoration examples 
in Figures 6(d) - (f), respectively. 

To quantify the stability of the JPDS estimates as 
a function of the number of realizations, we consider 
the sample variance about the sample mean. This 
is a measure of internal consistency: if the variance 
of restorations formed from disjoint sets of data de- 
creases as the number of realizations increases, then 
we conclude that restorations become more consistent 
with more realizations. If this variation is also quite 
small, then we can also conclude that any random 
false detail is too small to be considered problematic. 
A common method for measuring the variance of an 
estimator about its mean is via the coefficient of vari- 
ation (Frieden 1983), the square of which is defined 
at each pixel as 


L P(x) 


( 21 ) 


where ft is the ^th sample from a set of L (= 100/ J), 


J-realization restorations, and 
1 L ~ 

/w = z E^w ■ 

L t = l 


( 22 ) 


Figure 7 plots the spatial average of c^(x), denoted 
by c v , as a function of the number of realizations. By 
averaging over the pixels, we characterize the varia- 
tion about the mean at a pixel with a single number. 
The error bars on this plot represent the confidence 
in c„ to within one standard deviation. The stan- 
dard deviation of c v was computed as the square root 
of the variance of equation (21), which takes into ac- 
count the correlation of c v (x) from pixel to pixel. We 
note that Figure 7 displays a monotonic behavior from 
which we conclude that estimates become more con- 
sistent as realizations are added. Also, the greatest 
gains in consistency are made for small numbers of 
realizations, and c v is less than 1% for J = 10. Be- 
yond this point the gains are small, supporting our 
visual assessment of the sample restorations in Fig- 
ure 6. Figure 7 does not indicate how the mean im- 
ages change as a function of number of realizations. 
The RMS error between any pair of mean images has 
an average of about e — 0.3%, which reflects a high 
consistency for the expected restoration regardless of 
the number of realizations used. 

5.4. External Consistency of Phase-Diverse 
Speckle Restorations 

In the previous section we cited past evidence and 
presented new results that confirm the internal consis- 
tency of phase-diverse speckle object estimates. Ex- 
ternal validation of phase aberration estimates from 
conventional phase diversity has been demonstrated 
via the fixed aberrations of the SVST. These aber- 
rations were estimated over time as the SVST turret 
position changed and were shown to evolve accord- 
ing to theoretical predictions (Lofdahl & Scharmer 
1994b). Speckle imaging provides another indepen- 
dent way to validate the object estimates. For this 
work, the Fourier amplitudes have been reconstructed 
with the Labeyrie (1970) method, whereas the phases 
have been estimated with the Knox and Thompson 
(1974) technique. The Fourier amplitudes were cal- 
ibrated with a model of the Earth’s atmosphere by 
Korff (1973). Fried’s (1966) parameter, the only free 
parameter of the atmospheric model, was estimated 
with the spectral ratio technique (von der Liihe 1984) 
to be r 0 = 18.5 cm. A detailed description of the im- 
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plementation of the above-mentioned techniques has 
been given by von der Liihe (1993). 

Figure 8(a) contains the speckle reconstruction, 
which should be compared with the 100-realization 
phase-diverse speckle estimates from JPDS and PPDS 
in (b) and (c), respectively. Note that the overall 
features are identical in all 3 restorations, including 
subtle intensity variations along the edges of some 
granules. There is no evidence of artifactual details 
in either of the phase-diverse speckle restorations. In 
fact, the speckle restoration appears to be slightly 
smoother, with less high-frequency detail. To rein- 
force the important point that 100 image pairs are 
many more than required for phase-diverse speckle, 
we show examples of JPDS and PPDS 5-realization 
restorations in Figures 8(d) and (e), respectively. It 
is clear that excellent restorations can be obtained 
with a factor of 20 fewer realizations than used for 
(b) and (c). Because a reliable speckle reconstruc- 
tion like that in (a) could not be obtained with so few 
images, phase-diverse speckle can be viewed as an im- 
portant method for improving both the spatial and the 
temporal resolution of ground-based observations. 

The restorations in Figures 8(b) and (c) differ from 
the speckle restoration in (a) with an error of e = 2.0% 
and have an error of only e = 0.3%' between them. 
Some of the disparity is due to differences at high 
spatial frequencies, but we have found that small dif- 
ferences exist at low spatial frequencies as well. Fig- 
ure 9 plots the radially-averaged power spectra of 
the 3 restorations. Each spectrum was formed by 
first applying a Hanning window to the restoration, 
performing a 2-D Fourier transform, and then tak- 
ing the magnitude-squared. After appropriate scal- 
ing, each spectrum was integrated over annuli in the 
spatial-frequency domain. As expected, the JPDS 
and PPDS spectra are very consistent out to 80% 
of the diffraction-limited cutoff frequency, and differ- 
ences beyond this point can be attributed to differing 
regularization techniques. We note that the power in 
the speckle reconstruction begins to depart from the 
others at about half of the cutoff frequency, and that 
the speckle reconstruction appears to have approxi- 
mately 30% more power than the others up to this 
point. 

Defining contrast as the ratio of the standard de- 
viation to the mean, we find that the contrast of the 
speckle restoration is higher (12.6%) than either of the 
phase-diverse speckle estimates (11.0%). One plausi- 
ble explanation for this could be due to an insufficient 


representation of the phase aberrations with only the 
first 15 Zernike polynomials. It has been demon- 
strated in simulations that with a zonal (pixel-by- 
pixel) parameterization of the pupil, the full contrast 
of the object can be recovered (Seldin, Paxman, fe El- 
ste 1995). These same simulations consistently under- 
estimated the contrast when using only the first 15 
Zernike polynomials. On the other hand, the speckle 
reconstruction relies on a model for the atmospheric 
turbulence and does not account for fixed telescope 
aberrations. So the possibility of an underlying model 
mismatch remains, and the correctness of the restored 
contrast remains an open issue that is best studied 
with controlled experiments. Despite this significant 
difference in contrast, we conclude that the phase- 
diverse speckle estimates are consistent with and pro- 
vide finer detail than the reconstruction obtained with 
an accepted speckle-imaging technique. 

5.5. Evidence of Anisoplanatism 

Space- variant blur is encountered when objects ex- 
tend beyond the isoplanatic patch associated with 
the intervening atmosphere. The imaging model 
used here is space-invariant, and any anisoplanatism 
should degrade the quality of the restorations. Phase- 
diverse speckle object reconstructions are suscepti- 
ble to geometric distortions associated with aniso- 
planatism because of the underlying assumption of 
a time-invariant object across all realizations. Fur- 
thermore, intuition suggests that JPDS wavefront es- 
timates, unlike PPDS wavefront estimates, could also 
be hurt by these inter- realization geometric' distor- 
tions. The fact that both PPDS and JPDS succeed 
so well here is due in part to the selection of an 
appropriately-sized image patch, which was guided by 
previous analysis of this data sequence by Lofdahl and 
Scharmer (1994b). 

Despite the successful restorations obtained with 
these images, there is still evidence of anisoplanatism 
in this sequence. To characterize the geometric- dis- 
tortion, each image, dj ki in the sequence was resam- 
pled (destretched) to match the estimated image. g jk 
(eq. [20]), formed from the 100-realization JPDS es- 
timate in Figure 8(b). The destretching for the jtli 
realization is estimated by computing local correla- 
tions between dj k and gj k on a 3 x 3 segmentation 
of the images for both the conventional and diversity 
channels. The 9 shifts obtained from the 3x3 grid - >f 
correlations are averaged over the two channels, and a 
coordinate transformation that defines the destreteh- 
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ing function is formed by interpolating between the 
average shifts at the center of each segment. The 
degree of distortion, or image warping, can be sum- 
marized with a single number 

W = > (23) 

where xj and yd are the RMS shifts in the two direc- 
tions computed over all segments. Thus, large values 
of w reflect large deviations from a global shift and a 
higher degree of warping across the image. 

Figure 10(a) plots w (in units of pixels) versus the 
fidelity metric, R , for each of the 100 JPDS phase 
aberration estimates. There is a clear correlation 
(0.74) between R and the amount of distortion - re- 
alizations with the worst distortion yield the worst 
fidelity metric. In contrast, there is no apparent cor- 
relation between the RMS of the estimated aberra- 
tions and R , indicating that the strength of the tur- 
bulence is not the primary source of poor estimates. 
A new sequence of destretched images, djk , was cre- 
ated, and the fidelity metric was recomputed via equa- 
tion (19) using djk in place of djk- Note that the same 
phase aberration estimates were used to compute the 
“destretched” fidelity metric, and that the most im- 
provement in this metric occurs for the cases with 
large distortion. Remarkably, the improvement in the 
average fidelity metric after destretching is 3 times 
greater than the corresponding improvement when 
moving from conventional phase diversity to phase- 
diverse speckle. Figure 10(b) also suggests that aber- 
ration estimates are not hurt by the joint-estimation 
approach at this level of anisoplanatism. Although 
the improvement in the fidelity metric is only a frac- 
tion of a percent, we conclude that all phase-diversity 
techniques would benefit from the accommodation of 
anisoplanatic effects, particularly in the presence of 
stronger turbulence. An obvious approach would be 
to apply JPDS or PPDS to destretched data. Al- 
ternatively, parameters in an anisoplanatic-imaging 
model could be directly estimated from the original 
data (Paxman, Thelen, & Seldin 1994). 

6. Conclusions 

In this investigation, we have evaluated phase- 
diversity techniques, including two implementations 
of these techniques, and provided credibility for the 
scientific utility of these relatively novel observational 
techniques for solar astronomy. Using data collected 
with only modest instrumentation, we have recon- 


structed near diffraction-limited images of solar gran- 
ulation. 

6.1. Comparison of Implementations 

The ERIM and SVST implementations of phase- 
diversity techniques have been compared. In spite of 
significant differences in noise model, estimator, op- 
timization algorithm, regularization, and edge treat- 
ment, the phase-diversity estimates were found to be 
virtually identical. The invariance of estimates un- 
der differing noise models can be understood by con- 
sidering the similarities of the Poisson and additive 
Gaussian noise models for this application. Poisson 
noise is proportional to the square root of the inten- 
sity at each detector element. In addition, when the 
intensity is sufficiently large, the cumulative distribu- 
tion function for each Poisson random variable is well- 
approximated by a normal distribution (Feller 1968). 
For low-contrast images such as aberrated and/or de- 
focused images of solar granulation, the Poisson noise 
will be approximately constant across the image, thus 
resembling additive Gaussian noise. In fact, our re- 
constructions suggest that the Gaussian noise mode! 
is entirely adequate for imaging solar granulation. 
However, the Poisson model is the more general of the 
two for describing photo-detection events, and it can 
also be expanded to accommodate additive Gaussian 
noise sources, such as CCD readout noise. A Poisson 
model may be important when using low-noise cam- 
eras to image high-contrast scenes such as sunspot 
umbrae with umbral dots or for photon-limited imag- 
ing scenarios such as narrow-band solar imaging or 
nighttime astronomy. 

An advantage of the Gaussian noise model is that 
it leads to an optimization search within a reduced- 
dimension parameter space. The resulting estimation 
algorithm, as currently implemented by the SVST 
group, requires significantly fewer operations than the 
current ERIM algorithm, when operating with a suf- 
ficiently small number of aberration parameters, .U 
The SVST algorithm computations are dominated by 
Fourier transform calculations. The total number of 
FFTs required to perform a phase-diversity recon- 
struction is / • (7 4- 4 M), where / is the number of it- 
erations. For the data analyzed here, iterations were 
stopped when the RMS of the change of the wave- 
front was less than 10“ 3 radians, which is quite con- 
servative. This stopping criterion led to an average of 
about / = 12. Quality reconstructions are obtained 
if iterations stop when the metric changes less than 
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0.5%, giving (/> = 4.7. With M - 12 and / = 5, it 
takes 275 J FFTs to process J realizations with the 
SVST program. 

The ERIM algorithm is also dominated by two- 
dimensional FFTs. Approximately (3V + 4) J K FFTs 
are computed per iteration, where V is the number 
of likelihood evaluations made during the conjugate- 
gradient line search. Because this expression is inde- 
pendent of the number of aberration parameters, M , 
there is no additional computational cost per iteration 
for a fine parameterization of the wavefront. For the 
data analyzed herein, V = 10 and typically the num- 
ber of iterations performed was I = 60. These pro- 
cessing criteria are very conservative. We have since 
shown that visually indistinguishable object estimates 
are produced when relaxing the line search to V = 5 
and stopping the iterations when the metric changes 
by 0.02%, giving approximately / — 40. Using V = 5, 
I = 40, and K — 2 w'e see that the ERIM code re- 
quires approximately 1520J FFTs, considerably more 
than is required for the SVST code. In cases with M 
sufficiently large, such as a zonal aberration param- 
eterization (i.e. using pixels in the pupil), the ERIM 
code probably has a computational advantage, owing 
to the independence of the FFT count per iteration 
on M. It should be remembered that neither of the 
algorithms has been truly optimized with respect to 
computational efficiency. 

Another implementation issue is whether estimates 
are made in a joint or a partitioned fashion. Object 
estimates were found to be virtually identical using 
JPDS and PPDS. However, JPDS aberration esti- 
mates were found to be slightly more consistent with 
the data than were the PPDS estimates. Although 
the small improvement in aberration estimates pro- 
duced by JPDS is of little practical value in the re- 
constructions shown, it does suggest that a joint- 
estimation strategy could be of value when imaging 
in the presence of stronger turbulence or in regimes 
of reduced signal strength. We also believe that im- 
provement due to joint estimation is currently lim- 
ited by anisoplanatic effects that vary from realiza- 
tion to realization and by the under-parameterization 
of wavefronts. Both ERIM (Seldin, Paxman, k El- 
ste 1995) and SVST groups have found evidence that 
15 Zernikes is an under-parameterization of the wave- 
fronts, but the optimal number of wavefront parame- 
ters remains an open issue. 

When comparing the two approaches to treating 
edge effects, we see that both methods effectively han- 


dle problems associated with objects that extend be- 
yond the FOV. The ERIM guard-band method allows 
for the reliable retrieval of object pixels up to the edge 
of the detector-limited FOV, providing more efficient 
use of camera pixels and rendering a method for utiliz- 
ing interesting phenomena recorded along the camera 
borders. 

6.2. Comparison of Techniques 

Although estimated wavefronts were found to be 
essentially the same for all phase-diversity techniques, 
our results reiterate that phase-diverse speckle object 
estimates are significantly better than phase-diversity 
estimates. This can be understood by considering the 
OTFs involved. The OTFs for any one realization 
tend to reach very low values at isolated spatial fre- 
quencies, and the inversion of single-realization data 
can give false detail at these spatial frequencies. How- 
ever, these troublesome spatial frequencies vary from 
realization to realization so that when multiple real- 
izations are used, the likelihood of restoring false de- 
tail is dramatically diminished. Furthermore, we have 
quantitatively shown that estimates become increas- 
ingly more consistent as the number of realizations is 
increased. The added value of additional realizations 
can only increase when anisoplanatism is accommo- 
dated. 

Conventional speckle imaging can also be com- 
pared with phase-diverse speckle. Speckle imaging 
has the hardware advantage that only a single speckle 
camera is required. On the other hand, simple op- 
tical designs have been implemented that put both 
diversity images on a single camera, at the expense 
of reduced FOV. Speckle imaging also offers a com- 
putational advantage, requiring only J FFTs of size 
V 2 . Phase-diverse speckle has an important advan- 
tage in that excellent object reconstructions can be 
derived from a relatively few (J « 5) realizations, 
whereas speckle imaging requires a much larger num- 
ber ( J « 100) of realizations. With fewer realizations 
required, phase-diverse speckle can more readily be 
used for time-series observation of the evolution of 
solar phenomena. Whereas conventional speckle re- 
quires a separate speckle calibration step that relies 
on an atmospheric model (von der Liihe 1993) that 
normally doesn’t accommodate fixed telescope aber- 
rations, phase-diverse speckle requires no such step 
and is not hurt by fixed aberrations. In fact, because 
phase-diverse speckle estimates individual phase aber- 
ration realizations, these estimates can be averaged 
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to determine fixed aberrations (Lofdahl & Scharmer 
1994b). 

6.3. Future Directions 

Although there are issues which deserve greater 
study, including optimal wavefront parameterization, 
the efFects of anisoplanatism, and the effects of atmo- 
spheric evolution during an exposure time, we con- 
clude that phase-diversity techniques are sufficiently 
mature to be used routinely for solar observations. A 
new phase-diversity beam-splitter that puts two di- 
versity channels on a single large CCD, giving each a 
useful FOV of approximately 700 x 1000 pixels, has 
recently been installed at the SVST. This is avail- 
able as a common-user instrument. Excellent data, 
consisting of time sequences (several hours in dura- 
tion) of active regions, have already been obtained. 
In addition, scientific data, consisting of images of 
bright points, pores, sunspots, plages, and network 
magnetic fields in white-light, Ha, and Mglbo that 
cover up to several hours have been recorded with a 
multi-camera system at the Vacuum Tower Telescope 
at Sacramento Peak. Preliminary analysis of both 
data sets confirm our expectation that phase-diverse 
speckle opens a new observational window, provid- 
ing excellent spatial and temporal resolution over ex- 
tended periods of time. 

We are also in the process of migrating phase- 
diverse speckle methods to faint-object regimes, in- 
cluding narrow-band observations in solar astronomy 
and nighttime astronomical observations. As an ex- 
ample, we have recently used phase-diverse speckle 
to successfully correct for residual aberrations in an 
adaptive-optics system when imaging binary stars 
(Seldin, Paxman, & Ellerbroek 1995). 
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Fig. 1. — Data-collection scheme for phase-diverse 
speckle. A conventional short-exposure image and 
an additional short-exposure image that is intention- 
ally defocused are simultaneously collected for each 
of multiple atmospheric realizations. 

Fig. 2.“ Examples of preprocessed data. Five real- 
izations of conventional (upper) and diversity (lower) 
data are shown. Examples of good seeing are shown 
in (a) and (b). The image pair in (c) represents an av- 
erage case, and (d) and (e) contain images from which 
conventional phase-diversity restorations are poor. 

Fig. 3. — Examples of phase-diversity object es- 

timates from ERIM (upper) and the SVST group 
(lower). The five restorations in (a) - (e) are derived 
from the corresponding data in Figure 2. 

Fig. 4. — Scatter plots of Zernike coefficients of esti- 
mated aberrations using 2 phase-diversity implemen- 
tations. The SVST and ERIM estimates for the entire 
sequence of 100 images for coefficients 4-15 are shown. 
The high correlation is strong evidence of implemen- 
tation invariance. 

Fig. 5. — Scatter plot (over realization) of the fidelity 
metric R(<f>j) for SVST-group phase diversity versus 
ERIM phase-diverse speckle. R reflects the accuracy 
of the aberration estimates, and we observe that the 
joint phase-diverse speckle approach yields a better 
fidelity metric in 97 out of 100 cases. 

Fig. 6. — Representative object restorations for dif- 
ferent numbers of realizations. The restorations in 
(a) - (f) are derived using J image pairs for J = 

I, 2, 5, 10, 25, 50, respectively. Note that we show only 
one of many candidate restorations for each value of 

J, but the same input image-pairs are always retained 
as J is increased. The restorations are quite consis- 
tent for J > 5. 

Fig. 7. — Spatially averaged coefficient of variation, 
c V) in the object estimate as a function of number 
of realizations, J. The error bars were computed us- 
ing the variance of c v and are one standard deviation 
above and below the estimated value. 


Fig. 8. — Speckle and phase-diverse speckle restora- 
tions. (a) 100-realization speckle-imaging restora- 
tion; (b) 100-realization ERIM joint phase-diverse 
speckle restoration; (c) 100-realization SVST-group 
partitioned phase diverse speckle restoration; (d) 5- 
realization ERIM restoration; (e) 5-realization SVST- 
group restoration. 


Fig. 9. — Radially-averaged power spectra for 

speckle, ERIM joint phase-diverse speckle (JPDS), 
and SVST-group partitioned phase-diverse speckle 
(PPDS) restorations. The diffraction-limited cutoff 
frequency has been normalized to unity. 


Fig. 10. — Scatter plots of degree of distortion, u . 
versus fidelity metric, K, before and after dest retch- 
ing of the data, (a) The fidelity metric for the HK) 
ERIM phase-diverse speckle aberration estimates is 
correlated with the distortion; (b) After destretching 
the images, these same estimates yield fidelity met- 
ric values that are smaller on average and much less 
correlated with the original distortion. 
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